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We describe the general aspects of Monte Carlo Collision Generators suitable for 
cosmic ray nucleon-Air and nuclei-Air interactions, including accelerator and col- 
lider data. The problem of the extrapolation at 3 energy decades above the LHC 
of the main features of high energy collisions is discussed and under theoretical 
and phenomenological assumptions, the properties of the longitudinal and lateral 
development of giant extensive air showers simulated with the CORSIKA program 
are presented. The determination of the primary energy near 10^'^ eV is examined 
for different observables, total size, densities of charged particles interpolated at 
600 m from shower core. 

The extensive air shower data collected around LHC energy is in better agreement 
with models of large multiplicities. Beyond this energy, the extrapolation carried 
assuming the diquark breaking mechanism can change the classic conversion to 
primary energy and such circumstance can have consequences on the validity of 
the GZK cut off. 

In those conditions, we have simulated large and giant air showers taking into 
account, in addition, new processes, such as diquark breaking, and topological 
problems involving adequate structure functions for lateral distributions, up to 
energies exceeding 10^" eV for P.AUGER and EUSO experiments. 



1 Monte Carlo Collision generators for Cosmic Rays and 
extrapolation at ultra high energy 

In order to simulate in a reasonable time the 4-dimensional development of Ex- 
tensive Air Showers at very high energy, it is important to elaborate collision 
generators reproducing rapidly the detailed features of multiple production 
observed in accelerators and colliders. Among the variety of models imple- 
mented in CORSIKA, there are microscopic and phenomenological models. 
The first ones include all the steps of the parton momenta generation from 
the parton distribution functions related to valence quarks and diquarks, sea 
quarks and gluons, convoluted with the fragmentation of the respective strings 
into secondary hadrons. The phenomenological models provide fast and di- 
rect hadron sampling, but take into account some global characteristics of the 



Gribov-Regge theories with parameters adjusted to the experimental data. 
One common feature is the separation between the non single diffractive com- 
ponent and the diffractive component (single and double) with the inclusion 
of soft and hard mechanisms. 

As an example (fig.|l|), the reproduction of the inclusive data with the model 
HDPM2 (2nd version of the hybrid dual parton model) with parameters tuned 
to fit the experimental data of FermiLab at ^/s = 630GeV gives a total av- 
erage inelasticity for p-p collisions of 0.7 instead of 0.5 when adjusted to the 
previous measurements of UA5. The forward trajectory on Fig. |^ is not con- 
strained above 5.5 units of pseudo rapidity and several models taken in option 
for CORSIKA can have different trajectories in the forward region with a good 
agreement to the observed data, but rather different inelasticities. A cosmic 
ray cascade is initiated by the unique interaction of the primary particle and 
the fluctuations of the elementary act have to be reproduced carefully taking 
into account the semi-inclusive data as well as the correlation between charged 
and neutral secondaries. The extrapolation at ultra high energy can be car- 
ried with the model HDPM2 taking into account recent features of collider 
physics such as pt versua.central rapidity density (UAl-MIMI exp.) and re- 
cent results of FermiLab El for pseudo-rapidity up to 5.5. The pseudo-rapidity 
distributions obtained □ with HDPM2 for 2000 collisions NSD are shown on 
Fig. 0. The distributions on Fig. |l| (histograms for HDPM2, fuU line for QS- 
JET model) give for HDPM2 different extrapolations at 10^" GeV following 
the PDF assumed {BO and B_). 

The correlation between the central rapidity density and the average trans- 
verse momentum < pt > plays also one role in the muon and hadron radial 
distributions: the dependence of the transverse momentum on energy turns 
to a permanent increase on Fig. |^. 

The p-p collision is transformed to p-A collision following Glauber's con- 
siderations and Nuclei- Air collisions are treated by the abrasion evaporation 
procedure. 

2 Simulation of giant EAS 

The electromagnetic cascade is treated by EGS4 involving cross sections, tak- 
ing into acount the LPM effect, for bremsstrahlung and pair production. The 
nuclear photo production has been implemented in the simulation and its 
contribution is compared versus the primary photon energy to the previous 
processes on the Fig. |^ □. 

The fluctuations of multiplicity are governed by the negative binomial 
distribution and at the average charged multiplicity is close to 
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Figure 1. Left: pseudo-rapidity distribution with HDPM2. 
Right: average pt vs. energy dependence with HDPM2. 



300 for HDPM (1000 for QGSJET); this corresponds to energy densities of 
30 GeV/fm^ in average... and we have not considered here the important 
modifications from a possible phase transition to QGP. 

The longitudinal developments obtained for proton, photon and iron induced 
extensive air showers are compared on Fig. We notice that the showers 
(here averaged on 10 individual cascades simulated with CORSIKA) have 
their maximum near 1000-1500 m altitude for proton and iron. The 7 shower 
initiated at 50 g/crri^ exhibits a different behaviour due to the LPM effect; a 
rule of thumb for hadronic showers is to divide by 2 the total energy in GeV 
to get Njnax- This maximum depth corresponds to a minimum of fluctua- 
tions, suggested by analytic cascade theory and confirmed by the Monte Carlo 
simulation, indicating the well adapted localisation of AGASA and AUGER 
experiments. The depth of the maximum depends on the logarithm of the 
primary energy, when N„iax remains proportional to energy with a similar 
factor for a large variety of interaction models. 

A more delicate problem is the estimation of Nmax] it can be read directly 
from the cascade curve derived from fluorescence measurements, as performed 




Figure 2. Left: Photo-production cross section in air compared with bremsstrahlung and 

pair production cross sections reduced by the LPM effect. 

Right: longitudinal development of EAS for 7, p and Fe primaries. 



by the Fly's Eye or Hires, but the approach with a small number of detectors 
hit at large distances (1 to 1.5 km) turns to a hopeless topological problem. 

3 Topological aspects of radial distribution 

The profile of the lateral distribution assumed and the method of core local- 
ization are especially important. We propose instead of NKG and other Euler 
Beta functions, the employment of the gaussian hypergeometric formalism 
giving also normalization and better skewness, under the form: 

f{x) = g(s)a;"-"(x + + dx)-'= (1) 

which has the advantage (for values of parameters respecting the conditions 
of convergence s — a + 2 > and c—2s + b — 2 > 0) to be exactly normalized in 
terms of Gaussian Hypergeometric function i^/fc = F(c, s—a+2, c+b—s; 1 — d) 
by: 

/ X ^ r{c + b-s) 

2nr{s-a + 2)r{c-2s + b + a-2)FHG 

The empirical distributions, such as AGASA function i, as underlined 
Vishwanath □ enters in the category of Hypergeometric Gaussian functions 



Table 1. Best parameters to simulated e+e + muons (all charged) lateral distribution fit 
using JNCOl formula. 
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under the general form of structure function 

f{x) ^ Ce ■ x-"^ ■ (1 + a;)-'"-") • (1 + dxY<^ (3) 

with the conditions 2 — a > and /3 + 77 — 2 > 0. The value used in AG AS A 
function for the coefRcient Ce is just an approximation; the exact value is 

^ _ r(/3 + ,7-a) 1 



27r • r(2 - a) • r(/3 + 77 - 2) Fhg 

where F^q = Fhg{Pi 2 — a, (3 + t] ~ a;l — d). The Hypergeometric Gaussian 
function can be easily calculated from the hypergeometric series; 



F^oia, b, c;z)^Y. c ^ 0, -1, -2, 

^ (c)„n! 

(a)„ = r(a + n)/r(a), (a)o = 1 



This equation is equivalent to our version (3) containing the age parameter s 
with the relations between respective coefficients: x = —, d = —, 
s = 1.03, a — a ~ s, rj = b — s + a (the value of s is taken from the 
longitudinal development simulated). 

We have adjusted with MINUIT the parameters of our hypergeometric func- 
tion (Table |l|) to the average lateral distributions (set JNCOl for charged, 
JNC02 for electrons) of groups of 10 showers at 10^° eV simulated with COR- 
SIKA (QGSJET model), as shown on Fig. |. 

In each case, the adjustment has been performed with 50 points from the 
simulation distributed from 0.1 m up to 10 km from axis position for charged 
particles (muons and electrons). The advantages of JNCOl formula can be 
seen on Fig. || and on Table |[ 




Figure 3. Left: fits to all charged particles lateral distribution from simulations (average 
from 10 EAS). Primary particle is a proton with energy lO'^^ GeV. Lines are normalized to 
^(600 m). 

Right: extra component to GZK's cut off in case of diquark breaJjing. Pull line: generated 
spectrum, dashed: reconstructed, versus logio{Eo/l GeV) 



The major part of the particles is contained inside 200 m from the axis 
and only the skewness of the hypergeometric function allows a reliable re- 
lation between size and density at 600 m. This function has been applied 
to localizations of the showers contained in the catalogues of Volcano Ranch 
and Yakutsk. The core position has been obtained by minimization with Mi- 
nuit program between different formulas available for lateral densities written 
versus the coordinates X, Y an 

Q{r) = gi^iX-X,r + {Y-Y,)^) (4) 

where the core coordinates Xc and Yc are taken as two additive parameters in 
the minimization. The adjustments are generally improved when compared to 
the original treatments, turning to lower sizes (in the case of Yakutsk formula) 



Table 2. Columns a present total number of charged particles A^e in 10^ , columns b 
the ratios Eo/Ne in GeV {Eo=l(}^'^ GeV) and columns c the ratios Q{600)/Ne in 10~* 
particles/m^. m(600) = at 600 m. 
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Table 3. Localization of the most energetic AGASA event: x^i £?(600), A^e and relative 
core distance from the original localization. Age parameter s fitted for A and B functions 
as 1.05 and 0.98, respectively. 
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and better approximation o£ the density at 600 m. The situation of the most 
energetic event of AGASA El is given in Table |. 

4 Fluorescence, Cerenkov and Radio emissions 

From the longitudinal development, we have developed a method of fast sim- 
ulation based on the structural stability of the subshowers (in the sense of 
catastroph theory). This is an efhcient alternative to the thinning technique. 
The lateral extension of the total amount of light received at a fixed distance 
from the axis is compared for both Cerenkov and fluorescence component and 



PRIMARY PROTON 1.e11 GeV 



Longitudinal Profile PRIMARY PROTON 3.e8 GeV 



■ 
■ 




square = C 
circle - 


;erenkov light 
=luorescence light 


• 

1 


















• 

• 
















■ 




4 


> 














■ 










■ 








■ 











2.5 5 7.5 10 12.5 15 17.5 20 

DISTANCE ( KM ) 













22 deg 


rees * 




• 








• 1 


■ ■ ■ 


• 




< 




■ 














■ 






■ 
















• 
■ 










square = 


average 


of 100 showers with DPM option = 


circle = average of 10 showers with 


diquark 


breaking 













200 400 600 800 1000 1200 
DEPTH ( g/cm^ ) 



Figure 4. 1. Cerenkov and fluorescence light for a proton at 10^'^ eV. - 2. Simulation with 
DPM model with option D (100 showers), and diquark breaking (10 showers). 



shown on Fig. 

The opportunity of CORSIKA which separates positive and negative parti- 
cles taught us that a regular negative excess close to 25% is present in the 
e.m. component of EAS, allowing the calculation of radio emission, following 
Askarian's effect. In Yakutsk experiment, it is possible to estimate the primary 
energy from the density at 600 m as well as from the Cerenkov component; in 
AUGER array, the direct calibration of showers detected simultaneously by 
the ground array and the fluorescence detector looks promising. 



5 Diquark breaking mechanism and GZK cut ofT 

The diquark breaking mechanism disturbs strongly the leading particle effect 
present in the different models used in cosmic rays. In the classical form of 
the dual parton model, the 3 valence quarks of the proton projectile are sepa- 
rated in a fast diquark and another valence quark slowed down. The diquark 
is recombined with one quark of the sea to produce, the most commonly, an 
outgoing leader baryon propagating the energy deeper in the cascade. The 3 
valence quarks separated will be recombined in various meson structures in 
pairs \ud >, \du >,... or neutral mesons as 1/V2{dd - uu). The confi gura- 
tion with the simultaneous final state for the valence quarks of 37r'^'s could 



be especially interesting with a probability of emergence that we can evaluate 
from the quark content and the quark additive model as 1/27. Such con- 
figuration (with intermediate final states of higher probabilities, one pair of 
charged pions and one neutral, one pair of neutral and one charged...) will 
transfer a large part of energy to the electromagnetic component and this 
energy will be definitely missing for both hadronic and penetrating compo- 
nents. Remembering that for the same primary energy, the cascade theory 
shows that one primary photon produces at maximum, approximately, two 
times more electrons, we can expect a large electron excess for some cascades 
initiated with diquark breaking. The longitudinal development calculated for 
protons of the same energy of 3 • 10* GeV is compared to a classical devel- 
opment, here the model HDPM2 with D-option, on Fig. ^. The electron size 
at maximum is doubled in the assumption of diquark biEeaking and relatively 
rare recombination simultaneously in 3 neutral pions Ell. We have simulated 
10'* showers, between 10 EeV and 50 EeV, where the generation, following 
the primary energy differential spectrum, is cut by brute force. The showers 
simulated with axis distributed randomly are treated following the method 
of AGASA and the primary spectrum reconstructed is superimposed on the 
original one. It appears an extra component above 50 EeV, due to the show- 
ers including the diquark mechanism, considered previously, and such artefact 
could appear as a GZK violation. With AUGER array, such showers could be 
recognized very easily by an energy 3 or 4 times larger from the fluorescence 
detector than estimated by the water tank detectors. We observed also that 
such artificial component could appear from spacial topological situations at 
50 EeV; for instance, with an axis falling near one detector, all the others are 
hit by densities near the low density threshold and the answer very unstable 
of those surrounding detectors can turn to an important overestimation of the 
primary energy. 



6 Conclusions 

The hypergeometric approach gives a better accuracy in the interpolation of 
densities at 600-1000 m, a more reliable estimation of the shower size (when 
the axis is in the array), a better axis localization and finally a more correct 
constraint of the primary energy. The position of the axis is crucial as the 
density near 600 m varies as r^'^ and an uncertainty of 50 m on axis turns to 
an error of 40% at least on primary energy. Topological difficulties and also 
non standard aspect of multiple production needs probably more attention to 
be able to rule out definitely the cut off of GZK. 
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